New Data (helpdat: HELP study)

Today’s main data set comes from the Health Evaluation and Linkage to Primary Care trial, and is stored as HELPrct in the mosaicData package.

HELP was a clinical trial of adult inpatients recruited from a detoxification unit. Patients with no primary care physician were randomized to receive a multidisciplinary assessment and a brief motivational intervention or usual care, with the goal of linking them to primary medical care.

The helpdat data

We will look at 453 subjects with complete data today.

helpdat <- tibble(mosaicData::HELPrct) |>
    select(id, cesd, age, sex, subst = substance,
           mcs, pcs, pss_fr)

df_stats(~ cesd + age + mcs + pcs + pss_fr, data = helpdat) |>
  gt() |> fmt_number(min:sd, decimals = 2) |>
  tab_options(table.font.size = 20)

helpdat categorical variables

helpdat |> tabyl(sex, subst) |> 
    adorn_totals(where = c("row", "col")) |>
    adorn_percentages(denominator = "row") |>
    adorn_pct_formatting() |>
    adorn_ns(position = "front") |>
    adorn_title(placement = "combined") |>
  gt() |> tab_options(table.font.size = 20)

Our Outcome (CES-Depression score)

p1 <- ggplot(helpdat, aes(sample = cesd)) +
  geom_qq(col = '#440154') + geom_qq_line(col = "red") + 
  theme(aspect.ratio = 1) +
    labs(y = "Sorted CES-D Scores", 
         x = "Standard Normal Distribution")

p2 <- ggplot(helpdat, aes(x = cesd)) +
  geom_histogram(aes(y = stat(density)), 
                 bins = 10, fill = '#440154', col = '#FDE725') +
  stat_function(fun = dnorm, 
                args = list(mean = mean(helpdat$cesd), 
                            sd = sd(helpdat$cesd)),
                col = "red", lwd = 1.5) +
  labs(y = "Number of Subjects", x = "CES-D Score")

p3 <- ggplot(helpdat, aes(x = cesd, y = "")) +
  geom_violin(fill = '#440154') +
  geom_boxplot(width = 0.3, col = '#FDE725', notch = TRUE, 
               outlier.color = '#440154') +
  labs(x = "CES-D Score", y = "")

p1 + (p2 / p3 + plot_layout(heights = c(4,1))) +
  plot_annotation(title = "CES-D Depression Scores from helpdat data",
       subtitle = "Higher CES-D indicates more depressive symptoms",
       caption = "n = 453, no missing data")

Describing our outcome CES-D

describe(helpdat$cesd)

helpdat$cesd : CESD at baseline 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
     453        0       58    0.999    32.85    14.23     10.0     15.2 
     .25      .50      .75      .90      .95 
    25.0     34.0     41.0     49.0     52.4 

lowest :  1  3  4  5  6, highest: 55 56 57 58 60

• Info measures the variable’s information between 0 and 1: the higher the Info, the more continuous the variable is (the fewer ties there are.)
• Gmd = Gini’s mean difference, a robust measure of variation. If you randomly selected two of the 453 subjects many times, the mean difference in cesd would be 14.23 points.

We have some labels in our data

str(helpdat)

tibble [453 × 8] (S3: tbl_df/tbl/data.frame)
 $ id    : int [1:453] 1 2 3 4 5 6 7 8 9 10 ...
  ..- attr(*, "label")= chr "subject ID"
 $ cesd  : int [1:453] 49 30 39 15 39 6 52 32 50 46 ...
  ..- attr(*, "label")= chr "CESD at baseline"
 $ age   : int [1:453] 37 37 26 39 32 47 49 28 50 39 ...
  ..- attr(*, "label")= chr "age (years)"
 $ sex   : Factor w/ 2 levels "female","male": 2 2 2 1 2 1 1 2 1 2 ...
  ..- attr(*, "label")= chr "sex"
 $ subst : Factor w/ 3 levels "alcohol","cocaine",..: 2 1 3 3 2 2 2 1 1 3 ...
  ..- attr(*, "label")= chr "primary substance of abuse"
 $ mcs   : num [1:453] 25.11 26.67 6.76 43.97 21.68 ...
  ..- attr(*, "label")= chr "SF-36 Mental Component Score"
 $ pcs   : num [1:453] 58.4 36 74.8 61.9 37.3 ...
  ..- attr(*, "label")= chr "SF-36 Physical Component Score"
 $ pss_fr: int [1:453] 0 1 13 11 10 5 1 4 5 0 ...
  ..- attr(*, "label")= chr "perceived social support by friends"

Scatterplot Matrix (code)

temp <- helpdat |>
    select(age, mcs, pcs, pss_fr, sex, subst, cesd)

ggpairs(temp)  ## ggpairs from the GGally package

We place the outcome (cesd) last (result on next slide.)

Saving the Data Set

write_rds(helpdat, "c05/data/helpdat.Rds")

Scatterplot Matrix (result)

Fitting using ols

The ols function stands for ordinary least squares and comes from the rms package, by Frank Harrell and colleagues. Any model fit with lm can also be fit with ols.

• To predict var_y using var_x from the my_tibble data, we would use the following syntax:

dd <- datadist(my_tibble)
options(datadist = "dd")

model_name <- ols(var_y ~ var_x, data = my_tibble,
                  x = TRUE, y = TRUE)

This leaves a few questions…

What’s the datadist stuff doing?

Before fitting an ols model to data from my_tibble, use:

dd <- datadist(my_tibble)
options(datadist = "dd")

Run (the datadist code above) once before any models are fitted, storing the distribution summaries for all potential variables. Adjustment values are 0 for binary variables, the most frequent category (or optionally the first category level) for categorical (factor) variables, the middle level for ordered factor variables, and medians for continuous variables. (excerpt from datadist documentation)

Why use x = TRUE, y = TRUE?

Once we’ve set up the summaries with datadist, we fit a model:

model_name <- ols(var_y ~ var_x, data = my_tibble,
                  x = TRUE, y = TRUE)

• ols stores additional information beyond what lm does
• x = TRUE and y = TRUE save even more expanded information for building plots and summarizing fit.
• The defaults are x = FALSE, y = FALSE, but in 432, we’ll want them saved.

Using ols to fit a model

Let’s try to predict our outcome (cesd) using mcs and subst

• Start with setting up the datadist
• Then fit the model, including x = TRUE, y = TRUE

dd <- datadist(helpdat)
options(datadist = "dd")

mod1 <- ols(cesd ~ mcs + subst, data = helpdat,
                 x = TRUE, y = TRUE)

Contents of mod1?

mod1

Linear Regression Model

ols(formula = cesd ~ mcs + subst, data = helpdat, x = TRUE, y = TRUE)

                Model Likelihood    Discrimination    
                      Ratio Test           Indexes    
Obs     453    LR chi2    295.10    R2       0.479    
sigma9.0657    d.f.            3    R2 adj   0.475    
d.f.    449    Pr(> chi2) 0.0000    g        9.827    

Residuals

      Min        1Q    Median        3Q       Max 
-25.43696  -6.74592   0.09334   6.16212  24.24842 

              Coef    S.E.   t      Pr(>|t|)
Intercept     55.3026 1.2724  43.46 <0.0001 
mcs           -0.6570 0.0337 -19.48 <0.0001 
subst=cocaine -3.4440 1.0055  -3.43 0.0007  
subst=heroin  -1.7791 1.0681  -1.67 0.0965  

New elements in ols

For our mod1,

• Model Likelihood Ratio test output includes LR chi2 = 295.10, d.f. = 3, Pr(> chi2) = 0.0000

The log of the likelihood ratio, multiplied by -2, yields a test against a \(\chi^2\) distribution. Interpret this as a goodness-of-fit test that compares mod1 to a null model with only an intercept term. In ols this is similar to a global (ANOVA) F test.

New elements in ols

Under the \(R^2\) values, we have g = 9.827.

• This is the \(g\)-index, based on Gini’s mean difference. If you randomly selected two of the subjects in the model, the average difference in predicted cesd will be 9.827.
• This can be compared to the Gini’s mean difference for the original cesd values, from describe, which was Gmd = 14.23.

Validate summaries from an ols fit

• Can we validate summary statistics by resampling?

set.seed(432)
validate(mod1)

          index.orig training    test optimism index.corrected  n
R-square      0.4787   0.4874  0.4737   0.0137          0.4650 40
MSE          81.4606  79.7851 82.2361  -2.4510         83.9116 40
g             9.8272   9.9133  9.8038   0.1095          9.7177 40
Intercept     0.0000   0.0000  0.2793  -0.2793          0.2793 40
Slope         1.0000   1.0000  0.9894   0.0106          0.9894 40

• The data used to fit the model provide an over-optimistic view of the quality of fit.
• We’re interested here in assessing how well the model might work in new data, using a resampling approach.

Build a nomogram for the ols fit

plot(nomogram(mod1))

Using the nomogram for the mod1 fit

Predicted cesd if mcs = 35 and subst = heroin?

Actual Prediction for this subject…

• The predict function for our ols fit provides fitted values.

predict(mod1, newdata = tibble(mcs = 35, subst = "heroin"))

       1 
30.52766 

• The broom package can also support rms fits

augment(mod1, newdata = tibble(mcs = 35, subst = "heroin"))

# A tibble: 1 × 3
    mcs subst  .fitted
  <dbl> <chr>    <dbl>
1    35 heroin    30.5

Assessing the Calibration of mod1

We would like our model to be well-calibrated, in the following sense…

• Suppose our model assigns a predicted outcome of 6 to several subjects. If the model is well-calibrated, this means we expect the mean of those subjects’ actual outcomes to be very close to 6.

Assessing the Calibration of mod1

We’d like to look at the relationship between the observed cesd outcome and our predicted cesd from the model.

• The calibration plot we’ll create provides two estimates (with and without bias-correction) of the predicted vs. observed values of our outcome, and compares these to the ideal scenario (predicted = observed).
• The plot uses resampling validation to produce bias-corrected estimates and uses lowess smooths to connect across predicted values.
• Calibration plots require x = TRUE, y = TRUE in ols.

Calibration Plot for mod1

set.seed(43299); plot(calibrate(mod1))

n=453   Mean absolute error=0.128   Mean squared error=0.02428
0.9 Quantile of absolute error=0.267

Influential Points for mod1?

The dfbeta value for a particular subject and coefficient \(\beta\) is the change in the coefficient that happens when the subject is excluded from the model.

which.influence(mod1, cutoff = 0.2)

$Intercept
[1] "8"   "351" "405" "433"

$mcs
[1] "351" "402" "450"

$subst
[1] "351"

• These are the subjects that have absolute values of dfbetas that exceed the specified cutoff (default is 0.2 but it’s an arbitrary choice.)

Residuals vs. Fitted plot from ols

plot(resid(mod1) ~ fitted(mod1))

Fitting all Residual Plots for mod1

To fit more complete residual plots (and other things) we can fit the lm version of this same model…

mod1_lm <- lm(cesd ~ mcs + subst, data = helpdat)

par(mfrow = c(2,2)); plot(mod1_lm); par(mfrow = c(1,1))

• Plots are shown on the next two slides. While the subject in row 351 is more influential than most other points, it doesn’t reach the standard of a problematic Cook’s distance.

First Two mod1_lm Residual Plots

Second Two mod1_lm Residual Plots

Non-Linear Terms

In building a linear regression model, we’re most often going to be thinking about:

• for quantitative predictors, some curvature…
  • perhaps polynomial terms
  • but more often restricted cubic splines
• for any predictors, possible interactions
  • between categorical predictors
  • between categorical and quantitative predictors
  • between quantitative predictors

Polynomial Regression

A polynomial regression involves a polynomial in the variable x of degree D (linear combination of powers of x up to D.)

• Linear: \(y = \beta_0 + \beta_1 x\)
• Quadratic: \(y = \beta_0 + \beta_1 x + \beta_2 x^2\)
• Cubic: \(y = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3\)
• Quartic: \(y = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3 + \beta_4 x^4\)

An orthogonal polynomial sets up a model design matrix and then scales those columns so that each column is uncorrelated with the others.

Use pcs to predict cesd?

• Let’s look at both a linear fit and a loess smooth to see if they indicate meaningfully different things about the association between pcs and cesd

ggplot(helpdat, aes(x = pcs, y = cesd)) + 
    geom_point(size = 2) +
    geom_smooth(method = "loess", formula = y ~ x, 
                se = FALSE, col = "blue") +
    geom_smooth(method = "lm", formula = y ~ x,
                se = FALSE, col = "red") + 
    labs(title = "Linear and Loess Fits for cesd vs. pcs")

Polynomial regression with ols

dd <- datadist(helpdat)
options(datadist = "dd")

mod_B1 <- ols(cesd ~ pcs, 
              data = helpdat, x = TRUE, y = TRUE)
mod_B2 <- ols(cesd ~ pol(pcs, 2), 
              data = helpdat, x = TRUE, y = TRUE)
mod_B3 <- ols(cesd ~ pol(pcs, 3),
              data = helpdat, x = TRUE, y = TRUE)

• Note the use of pol() from the rms package here to fit orthogonal polynomials, rather than poly() which we used in an lm fit.

Model B1 (linear in pcs)

mod_B1

Linear Regression Model

ols(formula = cesd ~ pcs, data = helpdat, x = TRUE, y = TRUE)

                 Model Likelihood    Discrimination    
                       Ratio Test           Indexes    
Obs      453    LR chi2     40.57    R2       0.086    
sigma11.9796    d.f.            1    R2 adj   0.084    
d.f.     451    Pr(> chi2) 0.0000    g        4.177    

Residuals

     Min       1Q   Median       3Q      Max 
-28.4116  -7.8036   0.6846   8.7917  29.3281 

          Coef    S.E.   t     Pr(>|t|)
Intercept 49.1673 2.5728 19.11 <0.0001 
pcs       -0.3396 0.0522 -6.50 <0.0001 

Model B2 (quadratic poly. in pcs)

mod_B2

Linear Regression Model

ols(formula = cesd ~ pol(pcs, 2), data = helpdat, x = TRUE, y = TRUE)

                 Model Likelihood    Discrimination    
                       Ratio Test           Indexes    
Obs      453    LR chi2     40.68    R2       0.086    
sigma11.9915    d.f.            2    R2 adj   0.082    
d.f.     450    Pr(> chi2) 0.0000    g        4.199    

Residuals

    Min      1Q  Median      3Q     Max 
-28.387  -7.750   0.591   8.634  29.697 

          Coef    S.E.   t     Pr(>|t|)
Intercept 46.4007 8.7967  5.27 <0.0001 
pcs       -0.2136 0.3867 -0.55 0.5809  
pcs^2     -0.0014 0.0041 -0.33 0.7424  

Model B3 (cubic polynomial in pcs)

mod_B3

Linear Regression Model

ols(formula = cesd ~ pol(pcs, 3), data = helpdat, x = TRUE, y = TRUE)

                 Model Likelihood    Discrimination    
                       Ratio Test           Indexes    
Obs      453    LR chi2     48.70    R2       0.102    
sigma11.8991    d.f.            3    R2 adj   0.096    
d.f.     449    Pr(> chi2) 0.0000    g        4.556    

Residuals

     Min       1Q   Median       3Q      Max 
-27.5245  -8.2651   0.7988   8.9004  27.4480 

          Coef     S.E.    t     Pr(>|t|)
Intercept -13.4076 22.8605 -0.59 0.5578  
pcs         4.1323  1.5825  2.61 0.0093  
pcs^2      -0.1010  0.0354 -2.85 0.0046  
pcs^3       0.0007  0.0003  2.83 0.0049  

Store the polynomial fits

First, we need to store the values. Since broom doesn’t play well with ols fits, so I’ll just add the predictions as columns

cesd_fits <- helpdat |>
    mutate(fitB1 = predict(mod_B1),
           fitB2 = predict(mod_B2),
           fitB3 = predict(mod_B3))

Plot the polynomial fits

ggplot(cesd_fits, aes(x = pcs, y = cesd)) +
    geom_point() +
    geom_line(aes(x = pcs, y = fitB1),
              col = "blue", size = 1.25) +
    geom_line(aes(x = pcs, y = fitB2),
              col = "black", size = 1.25) +
    geom_line(aes(x = pcs, y = fitB3),
              col = "red", size = 1.25) +
    geom_text(x = 18, y = 47, label = "Linear Fit", 
              size = 5, col = "blue") +
    geom_text(x = 18, y = 39, label = "Quadratic Fit", 
              size = 5, col = "black") +
    geom_text(x = 18, y = 26, label = "Cubic Fit", 
              size = 5, col = "red") +
    labs(title = "Linear, Quadratic and Cubic Fits for cesd using pcs") 

Plot polynomial fits with Predict()

p1 <- ggplot(Predict(mod_B1)) + ggtitle("B1: Linear")
p2 <- ggplot(Predict(mod_B2)) + ggtitle("B2: Quadratic")
p3 <- ggplot(Predict(mod_B3)) + ggtitle("B3. Cubic")

p1 + p2 + p3

Types of Splines

• A linear spline is a continuous function formed by connecting points (called knots of the spline) by line segments.
• A restricted cubic spline is a way to build highly complicated curves into a regression equation in a fairly easily structured way.
• A restricted cubic spline is a series of polynomial functions joined together at the knots.
  • Such a spline gives us a way to flexibly account for non-linearity without over-fitting the model.

How complex should our spline be?

Restricted cubic splines can fit many different types of non-linearities. Specifying the number of knots is all you need to do in R to get a reasonable result from a restricted cubic spline.

The most common choices are 3, 4, or 5 knots.

• 3 Knots, 2 degrees of freedom, allows the curve to “bend” once.
• 4 Knots, 3 degrees of freedom, lets the curve “bend” twice.
• 5 Knots, 4 degrees of freedom, lets the curve “bend” three times.

Restricted Cubic Splines with ols

Let’s consider a restricted cubic spline model for cesd based on pcs with:

• 3 knots in modC3, 4 knots in modC4, and 5 knots in modC5

dd <- datadist(helpdat)
options(datadist = "dd")

mod_C3 <- ols(cesd ~ rcs(pcs, 3), 
              data = helpdat, x = TRUE, y = TRUE)
mod_C4 <- ols(cesd ~ rcs(pcs, 4), 
              data = helpdat, x = TRUE, y = TRUE)
mod_C5 <- ols(cesd ~ rcs(pcs, 5),
              data = helpdat, x = TRUE, y = TRUE)

Model C3 (3-knot spline in pcs)

mod_C3

Linear Regression Model

ols(formula = cesd ~ rcs(pcs, 3), data = helpdat, x = TRUE, y = TRUE)

                 Model Likelihood    Discrimination    
                       Ratio Test           Indexes    
Obs      453    LR chi2     40.79    R2       0.086    
sigma11.9901    d.f.            2    R2 adj   0.082    
d.f.     450    Pr(> chi2) 0.0000    g        4.206    

Residuals

     Min       1Q   Median       3Q      Max 
-28.3462  -7.7005   0.5098   8.6376  29.8454 

          Coef    S.E.   t     Pr(>|t|)
Intercept 47.3631 4.7053 10.07 <0.0001 
pcs       -0.2908 0.1187 -2.45 0.0146  
pcs'      -0.0624 0.1363 -0.46 0.6471  

Model C4 (4-knot spline in pcs)

mod_C4

Linear Regression Model

ols(formula = cesd ~ rcs(pcs, 4), data = helpdat, x = TRUE, y = TRUE)

                 Model Likelihood    Discrimination    
                       Ratio Test           Indexes    
Obs      453    LR chi2     51.31    R2       0.107    
sigma11.8648    d.f.            3    R2 adj   0.101    
d.f.     449    Pr(> chi2) 0.0000    g        4.590    

Residuals

     Min       1Q   Median       3Q      Max 
-28.3147  -8.2830   0.8559   8.8866  26.5458 

          Coef    S.E.   t     Pr(>|t|)
Intercept 33.3298 6.5742  5.07 <0.0001 
pcs        0.1464 0.1856  0.79 0.4308  
pcs'      -1.4383 0.4497 -3.20 0.0015  
pcs''      6.2561 1.9076  3.28 0.0011  

Model C5 (5-knot spline in pcs)

mod_C5

Linear Regression Model

ols(formula = cesd ~ rcs(pcs, 5), data = helpdat, x = TRUE, y = TRUE)

                 Model Likelihood    Discrimination    
                       Ratio Test           Indexes    
Obs      453    LR chi2     54.64    R2       0.114    
sigma11.8345    d.f.            4    R2 adj   0.106    
d.f.     448    Pr(> chi2) 0.0000    g        4.744    

Residuals

    Min      1Q  Median      3Q     Max 
-29.396  -7.928   1.016   8.762  26.974 

          Coef    S.E.   t     Pr(>|t|)
Intercept 39.0631 7.8282  4.99 <0.0001 
pcs       -0.0436 0.2332 -0.19 0.8517  
pcs'      -0.2952 1.0079 -0.29 0.7697  
pcs''     -3.1835 4.8079 -0.66 0.5082  
pcs'''    14.4216 8.3721  1.72 0.0857  

Plot all six fits?

p1 <- ggplot(Predict(mod_B1)) + ggtitle("B1: Linear")
p2 <- ggplot(Predict(mod_B2)) + ggtitle("B2: Quadratic")
p3 <- ggplot(Predict(mod_B3)) + ggtitle("B3. Cubic")
p4 <- ggplot(Predict(mod_C3)) + ggtitle("C3: 3-knot RCS")
p5 <- ggplot(Predict(mod_C4)) + ggtitle("C4. 4-knot RCS")
p6 <- ggplot(Predict(mod_C5)) + ggtitle("C5. 5-knot RCS")

(p1 + p2 + p3) / (p4 + p5 + p6)

Next Time

• Determining how to spend degrees of freedom on non-linearity
• The HERS data
• Fitting a more complex linear regression model
• Adding missing data into all of this, and using multiple imputation

How complex should our spline be?

Restricted cubic splines can fit many different types of non-linearities. Specifying the number of knots is all you need to do in R to get a reasonable result from a restricted cubic spline.

The most common choices are 3, 4, or 5 knots.

• 3 Knots, 2 degrees of freedom, allows the curve to “bend” once.
• 4 Knots, 3 degrees of freedom, lets the curve “bend” twice.
• 5 Knots, 4 degrees of freedom, lets the curve “bend” three times.

A simulated data set

set.seed(4322024)

sim_data <- tibble(
    x = runif(250, min = 10, max = 50),
    y = 3*(x-30) - 0.3*(x-30)^2 + 0.05*(x-30)^3 + 
        rnorm(250, mean = 500, sd = 70)
)

head(sim_data)

# A tibble: 6 × 2
      x     y
  <dbl> <dbl>
1  31.9 581. 
2  10.4 -90.5
3  38.0 476. 
4  15.6 232. 
5  39.6 474. 
6  42.5 721. 

The sim_data, plotted.

Fitting Non-Linear Terms with lm

We’ll fit:

• a linear model
• two models using orthogonal polynomials (poly()), and
• three models using restricted cubic splines (rcs())

sim_linear <- lm(y ~ x, data = sim_data)
sim_poly2  <- lm(y ~ poly(x, 2), data = sim_data)
sim_poly3  <- lm(y ~ poly(x, 3), data = sim_data)
sim_rcs3   <- lm(y ~ rcs(x, 3), data = sim_data)
sim_rcs4   <- lm(y ~ rcs(x, 4), data = sim_data)
sim_rcs5   <- lm(y ~ rcs(x, 5), data = sim_data)

augment() for fitted values and residuals

sim_linear_aug <- augment(sim_linear, sim_data)
sim_poly2_aug <- augment(sim_poly2, sim_data)
sim_poly3_aug <- augment(sim_poly3, sim_data)
sim_rcs3_aug <- augment(sim_rcs3, sim_data)
sim_rcs4_aug <- augment(sim_rcs4, sim_data)
sim_rcs5_aug <- augment(sim_rcs5, sim_data)

This will help us to plot the fits for each of these six models.

Add the Polynomial Fits

Restricted Cubic Spline Fits

Load and Partition Data

This is from the c4im example we used last Thursday.

c4im <- read_rds("c04/data/c4im.Rds")

c4im <- c4im |>
  mutate(fruit_c = fruit_day - mean(fruit_day))

set.seed(432)    ## for future replication
c4im_split <- initial_split(c4im, prop = 3/4)
train_c4im <- training(c4im_split)
test_c4im <- testing(c4im_split)

Fitting Restricted Cubic Splines with lm and rcs

For most applications, three to five knots strike a nice balance between complicating the model needlessly and fitting data pleasingly. Let’s consider a restricted cubic spline model for 1000/bmi based on fruit_c again, but now with:

• in temp3, 3 knots, and
• in temp4, 4 knots,

temp3 <- lm(1000/bmi ~ rcs(fruit_c, 3), data = train_c4im)
temp4 <- lm(1000/bmi ~ rcs(fruit_c, 4), data = train_c4im)

Spline models for bmi and fruit_c

Let’s try an RCS with 4 knots

m_4 <- lm(1000/bmi ~ rcs(fruit_c, 4) + exerany + health,
          data = train_c4im)

m_4int <- lm(1000/bmi ~ rcs(fruit_c, 4) + exerany * health,
          data = train_c4im)

m_4int coefficients

m_4int Residual Plots

par(mfrow = c(1,2)); plot(m_4int, which = c(1,2))

m_4int Residual Plots

par(mfrow = c(1,2)); plot(m_4int, which = c(3,5))

m_4 and m_4int in test sample

After back-transformation of fitted values of 1000/bmi to bmi:

bind_cols(testing_r2 |> select(model, rsquare = .estimate), 
          testing_rmse |> select(rmse = .estimate),
          testing_mae |> select(mae = .estimate)) |> 
  gt() |> fmt_number(columns = -model, decimals = 4) |>
  tab_options(table.font.size = 20)